Inertial range scaling in numerical turbulence with hyperviscosity 
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Numerical turbulence with hyperviscosity is studied and compared with direct simulations using 
ordinary viscosity and data from wind tunnel experiments. It is shown that the inertial range scal- 
ing is similar in all three cases. Furthermore, the bottleneck effect is approximately equally broad 
(about one order of magnitude) in these cases and only its height is increased in the hyperviscous 
case-presumably as a consequence of the steeper decent of the spectrum in the hyperviscous sub- 
range. The mean normalized dissipation rate is found to be in agreement with both wind tunnel 
experiments and direct simulations. The structure function exponents agree with the She-Leveque 
model. Decaying turbulence with hyperviscosity still gives the usual t~^''^^ decay law for the kinetic 
energy, and also the bottleneck effect is still present and about equally strong. 



I. INTRODUCTION 

In recent years there has been growing awareness of 
the detailed structure of the kinetic energy spectrum of 
hydrodynamic turbulence. In addition to the basic Kol- 
mogorov k~^/^ spectrum with an exponential dissipation 
range there are strong indications of intermittency cor- 
rections (possibly throughout the entire inertial range) 
and there is also the so-called bottleneck effect 

QH, i e. 

a shallower spectrum near the beginning of the dissipa- 
tive subrange; see also Ref. These features can be 
seen both in high resolution simulations and in mea- 
surements of wind tunnel turbulence . 

Over the past few years it has become evident that 
in numerical turbulence the bottleneck effect is rather 
pronounced ii, _7J. However, some of the simulations 
used hyperviscosity or other kinds of subgrid scale mod- 
eling. Hyperviscosity has frequently been used in turbu- 
lence studies in order to shorten the dissipative subrange 
0,0,0,0,01. However, hyperviscosity has also been 
suggested as a possible source of an artificially enhanced 
bottleneck effect Meanwhile, the apparent dis- 

crepancy in the strength of the bottleneck effect between 
simulations and experiments has been identified as being 
due to the difference in the diagnostics: in wind tunnel 
experiments one is only able to measure one-dimensional 
(longitudinal or transversal) energy spectra, while in sim- 
ulations one generally considers shell integrated three- 
dimensional spectra. The two are related by a simple 
integral transformation 0,0,0]. It turns out that, 
while the bottleneck effect can be much weaker or even 
completely absent in the one-dimensional spectrum, it is 
generally much stronger in the three-dimensional spec- 
trum [l||. 

In order to see the bottleneck effect in simulations, it is 



'Electronic address: nils.haugen@phy s.ntnu.no| 
t Electronic address: brandenb@nordita.dkJ 



important to have sufficiently large resolution of around 
1024'^ meshpoints. This raises the question to which ex- 
tent the bottleneck effect seen in simulations with hy- 
perviscosity is an artifact or a real feature that becomes 
noticeable only above a certain resolution. It is thus pos- 
sible that the reason for an exaggerated bottleneck effect 
in the hyperviscous simulation is related to the fact that 
hyperviscosity increases the effective resolution beyond 
the threshold above which the bottleneck effect can be 
seen. 

In this paper we consider forced hydrodynamic turbu- 
lence using hyperviscosity proportional to (instead of 
the usual viscosity operator). We find that the bottle- 
neck effect is enhanced in amplitude-but not in width, 
compared with direct simulations at the currently largest 
resolution of 4096'^ on the Earth Simulator One of 
the important results of these very high resolution simu- 
lations is that an inertial range begins to emerge that is 
clearly distinct from the bottleneck effect. Furthermore, 
the (negative) slope in the inertial range is steeper than 
the standard Kolmogorov power law exponent of 5/3 by 
about 0.1, so it is approximately 1.77. 

As in earlier papers , we consider weakly compress- 
ible turbulence using an isothermal equation of state. 
The root mean square Mach number is between 0.12 and 
0.13; for this type of weakly compressible simulations, we 
find that the energies of solenoidal and potential compo- 
nents of the flow have a ratio £^pot/£'soi ~ 10^'*-10^^ for 
most scales; only towards the Nyquist frequency the ratio 
increases to about 0.1. Compressibility is therefore not 
expected to play an important role. 



II. BASIC EQUATIONS 

We solve the compressible Navier-Stokes equations, 
Dn 1 



(1) 
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where D/Dt = d/dt + it • V is the advective derivative, 
p is pressure, p is the density, / is an isotropic, random, 
nonhehcal forcing function with power in a narrow band 
of wavenumbers, and 



P ^ 



(2) 



is the viscous force. Here, 

S(») = (_v2)"-is (3) 
is a higher order traceless rate of strain tensor, 

Sy = \ (ui^j + Uj,i) - 1% V • u (4) 

is the usual traceless rate of strain tensor, and commas 
denote partial differentiation. In the following we restrict 
ourselves to the case where fin = pt^n = const. Using the 
product rule, we can then rewrite Eq. |5| in the form 

Fvisc = (-1)"-!^ (v2"m + iv^f"-!) VV • w) . (5) 

For n = 1 we recover the normal diffusion operator for 
compressible flows. In the present paper we choose n — 3, 
so equation ||SJ) reduces to 

Fvi,e = ^(V6t^+|V4VV.M). (6) 

In the incompressible case, which is usually considered, 
the second term in Eq. vanishes. However, in the 
compressible case considered here this term is important 
to ensure momentum conservation. The local rate of ki- 
netic energy dissipation per unit mass is 



e = 2^3 (V^S)' 



(7) 



which is positive definite. 

We consider an isothermal gas with constant sound 
speed Cs, so that the pressure is given by p = c^p and 
p~^Vp = CgVlnp. The density obeys the continuity 
equation. 



Dlnp 
Dt 



= -V • u. 



(8) 



For all our simulations we have used the Pencil Code 
, which is a grid based high order code (sixth order in 
space and third order in time) for solving the compress- 
ible hydrodynamic equations. 



III. RESULTS 

We have calculated a series of models with resolutions 
varying between 64^ and 512'^ meshpoints using a third 
order hyperviscosity (n = 3). When changing the reso- 
lution, we keep the grid Reynolds number, here defined 
as 



rid 



(9) 



approximately constant. Here, ^Ny = 7r/5a; is the 
Nyquist wavenumber and 5x is the mesh spacing. Thus, 
when doubling the number of meshpoints, we can de- 
crease the viscosity by a factor of about 2^ = 32. This 
shows that hyperviscosity can allow a dramatic increase 
of the Reynolds number based on the scale of the box. 

Higher order hyperviscositics (rt = 3 and larger) have 
been studied previously 12, 20], but for us n = 3 is a 
practical limit, because we have restricted the maximum 
stencil length of all derivative schemes to three (in each 
direction), which is required for sixth order finite differ- 
ence schemes for our first and second derivatives fillj . 

In the following we consider the convergence of the 
energy spectrum for our hyperviscous simulations and 
compare with direct simulations. We discuss then the 
Reynolds number dependence of the normalized mean 
dissipation rate, and present finally the scaling behavior 
of the structure functions. Our basic conclusion is that 
in hyperviscous and direct simulations, as well as in wind 
tunnel experiments, the inertial range scaling is virtually 
identical and the width of the bottleneck is similar. 



A. Energy spectra 

Here and below we have calculated the energy dissipa- 
tion rate from the energy spectrum via 



e = 2vn / kl^E{k) dfc. 



(10) 



Here we have taken into account that in the code we 
employ a finite difference scheme which has always a 
discretization error, so we have to use the effective 
wavenumber in the expression above. The effective 
wavenumber is usually less than the actual one; see fig- 
ure 9.1 of Ref. |23|. For example, for the sixth order 
finite difference scheme, an analytic expression for k'^^ 
was given in Ref. p^ . while in the present case we have 



20 — 30 cos K + 12 cos 2k — 2 cos 3k, 



(11) 



where k = k5x is the wavenumber scaled by the mesh 
spacing 8x. Using the effective wavenumber becomes par- 
ticularly important in the hyperviscous case in order not 
to overestimate the contribution to e in Eq. H10() . 

The dissipation wavenumber, fed, is calculated from the 
relations e — k^u^^ and e = Vnk'^u\^. This leads to 



, 6ri-2 



kf forn 



3) 



(12) 



Again, for n = 1 one recovers the usual relation k^ — 
(tjv^YI^ . For larger values of n we find that, in order to 
make the location of the inertial range in direct and hy- 
perviscous simulations agree, we have to use an effective 
wavenumber /cd,eff that is larger than fed by a factor that 
is around 4 in our case, i.e. fed, off ~ 4fcd. 

In Fig.^we show the convergence of the energy spectra 
of hyperviscous runs for increasing resolution up to 512'^ 
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FIG. 1: Time- averaged energy spectra compensated by 
curves correspond to four different resolu- 
tions. All runs are with hyperviscosity. 
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FIG. 2: Time-averaged energy spectra, compensated by 
/j5/3^-2/3^ for the direct simulation with 4096^ meshpoints 
at ReA = 1201 (solid line) from figure 5 of Ref. Q and our 
hyperviscous simulation with 512'^ meshpoints (dashed line). 
Note that the bottleneck has a higher amplitude in the hy- 
perviscous case, but the inertial range has the same slope as 
for the simulation with 4096^ meshpoints. Our hyperviscous 
energy spectrum is scaled by a factor 1.1 in order to make it 
fall on top of the 4096'' result, i.e. our Kolmogorov constant 
is 1.1 times smaller than for the 4096^ simulation. 



meshpoints. All spectra are compensated by a k^/^^-"^/^ 
factor and the abscissa is normalized to the effective dis- 
sipation wavenumber A:d,eff- AH runs agree in the shape 
of the bottleneck and the subsequent dissipation sub- 
range, but the length of the inertial range varies from 
non-existent to about one order of magnitude. 

We now compare our 512'^ meshpoints hyperviscous 
run with the direct simulations of Kaneda et al. 0| on 
the Earth Simulator using 4096'^ meshpoints; see Fig. |21 




0.1 I ^^^^^ I 
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FIG. 3: Our simulation with 1024'^ meshpoints and normal 
viscosity show a bottleneck very similar to the bottlenecks in 
Ref , but due to lack of resolution we do not see any inertial 
range. 

We see that in both cases the bottleneck sets in at 
k/kd,cs ~ 0.03 and spans approximately one decade, but 
the dissipation subrange is longer in the direct simula- 
tions. The height of the bottleneck increases with in- 
creasing order of the hyperviscosity which is not 
surprising given that the steepness of the dissipative sub- 
range is the reason for the bottleneck effect in the first 
place 0. In agreement with Kaneda et al. 0|, we find 
that the slope of the energy spectrum in the inertial range 
is consistent with the fc"^ ''^ law found in the direct sim- 
ulation. The Kolmogorov constant is however slightly 
smaller (about xl.l) in our hyperviscous case. 

We should emphasize that, although we solve the com- 
pressible equations using finite differences, our direct 
simulations agree favorably with those using spectral 
methods solving the incompressible equations. This is 
shown in Fig. |3 where we compare simulations using 
1024^ meshpoints and normal viscosity with those of 
Ref H. These data have previously been discussed in 
R,efs. jld l2^ in connection with the bottleneck effect in 
hydrodynamics and hydromagnetic turbulence. 

We now compare with the data from a wind tun- 
nel experiment. Ideally we would like to translate the 
one-dimensional wind tunnel data into three-dimensional 
data |lS| | , but this involves differentiation which amplifies 
the noise in the data. Therefore we now compare one- 
dimensional energy spectra of our largest hyperviscous 
simulation with the energy spectrum from a wind tunnel 
experiment; see Fig. 01 We see that in our simulation the 
bottleneck has larger amplitude than in the wind tunnel 
experiment, but the (negative) slope of the inertial range 
spectrum is comparable in the two cases, i.e. 1.77. The 
Kolmogorov constant on the other hand is smaller by a 
factor of 1.5 in the hyperviscous case compared to the 
wind tunnel experiment. 

We feel that the value of a slope of 1.77 should be taken 
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FIG. 4: One-dimensional time-averaged energy spectra of our 
largest run with hyperviscosity compared with wind tunnel 
data with Re^ = 730 '5|. We have multiplied our energy 
spectra by 1.5 in order to make it fall on top of the wind 
tunnel data. 



FIG. 5: Plot of Ce as a function of Re^ for runs with third 
order hyperviscosity (n = 3). Triangles and plus signs rep- 
resent Reynolds numbers calculated based on Eq. 11611 with 
Rsao ~ 7.5 and 16, respectively, while for the plus signs 
Eq. I|14|l together with Eq. H15|l have been used. 



with caution, because it departs rather markedly from 
the value 1.70 expected from the She-Leveque relation 
p^ . Given that the inertial range is still relatively short, 
a slope of 1.70 can certainly not be excluded. 

It is customary to quote the Reynolds number based 
on the Taylor microscale j25| . 

A = \/5 Urms/t^rms- (13) 

Furthermore, Urms and Wrms are the rms velocity and 
vorticity, respectively. One usually takes the one- 
dimensional rms velocity for defining the Reynolds num- 
ber. 

Rex = uiuX/i^, (14) 

where — ^uf.^-^^. The wind tunnel experiments have 
ReA = 730. 

In the hyperviscous case the straightforward defini- 
tion of the Taylor microscale Reynolds number would 
be Rca = MiDA^/i^3, but this would lead to rather large 
values (~ 10^) which would not be meaningful in this 
context. Instead we define an effective viscosity from the 
actual mean dissipation rate and the modulus of the or- 
dinary rate of strain matrix, 

z^eff = (e)/(2S2), (15) 

which is then used to estimate the value of ly in Eq. (|14|l . 
In this way we find Rca = 340 for our largest simulation. 
Comparing with the high resolution direct simulations 
(Fig. O and with wind tunnel data (Fig. ^ we see that 
Rba = 340 probably is an underestimate for our hyper- 
viscous simulations. 

Alternatively one can define Rca as a measure of the 
width of the inertial range. Using relations that are valid 



in the standard case with 71 = 1, we have fcd,cff/^f ~ 
Re^/'' and Rca Re^/^ which yields 

ReA«RcAo(^)'^', (16) 

where we have introduced Rcao as a calibration parame- 
ter, and fcf is the forcing wavenumber or, more generally, 
the wavenumber of the energy carrying scale. If we set 
Rgao ~ 7.5, we can reproduce the result Rca = 340 for 
our largest run. On the other hand, if we choose to cali- 
brate Rcao such that our run with 512'^ meshpoints and 
the wind tunnel experiments have the same Rca = 730 
(see; Fig.QJ then we find Rcao — 16, which is perhaps a 
more reasonable estimate. 



B. Energy dissipation rate 

According to the Kolmogorov phenomenology, the 
spectral energy flux should be independent of k in the 
inertial range and equal to both the rate of energy input 
at large scales and the rate of energy dissipation at small 
scales. The constant of proportionality is of fundamental 
interest in turbulence research and one wants to know 
whether this value is independent of Reynolds number 
4, 25j. It is customary to define this coefficient as 

= {e)L/ul^, (17) 

where (e) is the mean energy dissipation rate and L the 
integral scale, which is usually defined as L = (37r/4)fc]~^, 
where 

fcfi== [ k-^E{k)dk [ E{k)dk, (18) 



5 



Structure functions 
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FIG. 6: Time-averaged total structure functions, [Sp ^ + 
2S'p'']/3, for p = 2 and p — 3. The two dotted liorizontal 
lines go through 0.7 and 1.0, confirming the expected scaling 
from the She-Leveque relationship. 
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FIG. 7: Structure function scaling exponents found using the 
concept of extended self similarity. We see that the longitudi- 
nal scaling exponents follow the She & Leveque scaling very 
well, while the transversal scaling exponents are somewhat 
more intermittent. 



is the spectrally weighted average of A: ^. 

The resulting normalized mean energy dissipation rate, 
calculated in this way, is shown in Fig. El as a function of 
Reynolds number. In the figure diamonds correspond to 
using Eq. p4|l and Eq. H15|l to find the Reynolds num- 
ber, while triangles and plus signs correspond to using 
Eq. ifTCjl with Rcao = 7.5 and 16, respectively. Figure 
shows that our results are in good agreement with both 
numerical 0, and experimental l Eq data. 



The spectral information can be supplemented by simi- 
lar scaling information in real space using structure func- 
tions. We define the longitudinal and transversal struc- 
ture functions 



S^^\r) = {{r ■ [u{x + r) - u{x)]y 



Sl'\r)^{{h-[u{x 



u{x)]y 



(19) 



(20) 



respectively. Here, r is the unit vector of r and fi is 
normal to r, so n • f = 0. The structure function of the 
three-dimensional velocity field, 



Spir)^{\u{x + r)-uix)\P), 
can then be written as 

5p(r) = i[5«(0 + 25W(r)]. 



(21) 



(22) 



We define the pth order structure function scaling expo- 
nent, C^p via the scaling relation 



Sp{r) (X r^p. 



(23) 



In Fig.Elwe plot the derivative of the double- logarithmic 
slope of the structure functions, dlnSp/dlnr. Inertial 
range scaling is indicated by a plateau in this graph. We 
find from the lower curve of Fig. that ^2 ~ 0.7. More 
importantly, in the upper curve of Fig. we show that 
53 (r) is consistent with linear scaling, i.e. ^3 = 1- Know- 
ing this we can use the extended self similarity to 
find the other structure function scaling exponents. The 
results are shown in Fig. [3 where we see that the lon- 
gitudinal structure function exponents follow the She & 
Leveque scaling very well, while the transversal struc- 
ture function is slightly more intermittent (i.e. the graph 
of versus p is more strongly bent). In particular we 
note that the extended self similarity gives C2 — 0.696 for 
the longitudinal component. 



D. Decaying turbulence 

Decaying isotropic turbulence is often considered an 
important benchmark of turbulence theories, and com- 
parisons with wind tunnel experiments and large eddy 
simulations are available We have carried out sim- 
ulations of decaying turbulence by stopping the driving 
at a time that will now be redefined to t = 0. In Fig. [S] 
we show that asymptotically 



(24) 



where uq is the initial [t — 0) rms velocity and r is a 
constant obtained from the intersection of the decay law 



extrapolated to {u 



It turns out that our law 
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FIG. 8: Mean squared velocity as a function of time for 
a simulation where the forcing has been stopped at t = to- 
In the 256"^ simulation we have Regrid ~ 0.20 (based on the 
initial value uq) and ruofci = 2.8 while in the 512^ simulation 
we have Regrid = 0.24 and ruoki = 3.1. The solid circles 
correspond to the experimental results of Kang et al. |2^. 



is consistent with n — 1.25, which is in agreement with 
recent wind tunnel data and large eddy simulations . 

The bottleneck effect is roughly unchanged; see Fig. El 
Its width is still about one order of magnitude in 
wavenumber, which is comparable to the experimental 
results li^. The height of the bottleneck is much less 
in the experimental data, because they show only a one- 
dimensional spectrum which gives a much weaker hump 
than the three-dimensional spectra p^ . In any case, we 
know already that the height of the bottleneck is artifi- 
cially enhanced by the use of hyperviscosity. 

In the present decay simulations one also sees the be- 
ginning of a subinertial range, leaving only a very short 
inertial range around 0.1 k(v-j,tY^^ ^ 0.3. The slope 
is compatible with the She-Leveque slope of 1.70, which 
corresponds to a residual slope of after compensat- 

ing with k^/^. Thus, there is no longer evidence for the 
more extreme correction of —0.1 suggested by the forced 
turbulence simulations H. 



IV. CONCLUSION 

The present investigations have shown that turbulence 
simulations with hyperviscosity are able to reproduce vir- 
tually the same inertial range scalings as simulations with 
ordinary viscosity. Specifically, the structure function ex- 
ponents show scaling behavior that is consistent with the 
She-Leveque model. However, the transversal struc- 
ture functions show a slightly higher degree of intermit- 
tency than the longitudinal ones. This, in turn, is quite 
consistent with a number of turbulence simulations by 



FIG. 9: Energy spectra for a decaying run. The abscissa 
is compensated by {usty^^ to make it dimensionless and to 
account for the slow decrease of the dissipation wavenumber. 
The ordinate is compensated by k^^'^ to show the location of 
the inertial range and by f^'''* to compensate for the decay. 



other groups l30l |. A possible explanation for the 
difference between longitudinal and transversal structure 
functions has been offered by Siefert & Peinke 0|, who 
find different cascade times for longitudinal and transver- 
sal spectra. The spectra show inertial range scaling sim- 
ilar to that found both in wind tunnel experiments Q 
and in very high resolution direct simulations In all 
three cases (hyperviscous and direct simulations as well 
as wind tunnel experiments) the inertial range spectrum 
is found to be compatible with the A:"^ *"^ behavior found 
by Kaneda et al. [Jl- As discussed above, this result is not 
compatible with the results from the structure function 
scalings and the She-Leveque relation. However, we be- 
lieve that the presently resolved inertial range is still too 
short to distinguish conclusively between 1.77 and the 
She-Leveque value of 1.70. Also, the simulation data of 
decaying turbulence suggest a weaker correction of 0.03, 
giving a slope of 1.70 that is compatible with the She- 
Leveque scaling. 

Another important result is that the width of the bot- 
tleneck seems to be independent of the use of hyperviscos- 
ity, and that only its height increases with the order of the 
hyperviscosity. This result is also confirmed in the case 
of decaying turbulence. Finally, we note that the nor- 
malized dissipation rate is independent of the Reynolds 
number, and that the asymptotic value of « 0.5 is in 
agreement with both experimental and numerical results 

One should of course always be concerned about the 
possible side effects of using hyperviscosity. One worry 
is that hyperviscosity ma y a ctually affect almost all of 
the inertial subrange Il4j. The current simulations 
confirm that the bottleneck requires at least an order of 
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magnitude in fc-space, and so does the dissipative sub- 
range, leaving almost no inertial range at all-even in a 
simulation with 1024^ meshpoints. Thus, using hypervis- 
cosity appears to be a reasonable procedure for gaining 
information about the inertial range at moderate cost, 
even though one should still use a reasonably high res- 
olution to isolate true inertial range features from those 
in the bottleneck subrange. On the other hand, hyper- 
viscosity is not a universally valid approximation. An 
example is in magnetohydrodynamics when magnetic he- 
licity is finite and a large scale magnetic field builds up 
in a closed or fully periodic box ;;33] . As long as it is pos- 
sible to understand the origin of peculiar features arising 
from hyperviscosity or hyper-resistivity (as is the case 
in helical hydromagnetic turbulence) there may well be 
circumstances where turbulence with hyperviscosity can 
provide a useful model for certain studies. One should 
bear in mind, however, that the height of the bottleneck 



depends on the order of the hyperviscosity. For example, 
if we choose n > 3 in Eq. (jSJ, the height of the bottleneck 
will be even more exaggerated [T^ . 
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